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RELAXATION  TECHNIQUES  FOR  THREE-DIMENSIONAL  STEADY  TRANSONIC  FLOW 
ABOUT  WINGS 

by  Li  Yimin 
ABSTRACT 

Bused  on  difference  equations  obtained  from  -small  perturbation  theory  and 
through  a  hyperbolic  tangent  transformation  which  maps  the  physieal  space  into  a  cube, 
the  three  dimensional  steady  transonic  flow  about  a  wing  is  computed.  Relaxation 
method  reduces  ihc  demand  on  computer  ni  -mor.v,  and  experience  is  gained  through  test 
runs.  Calculated  results  are  in  good  agreement  with  wind  tunnel  testa. 

In  transonic  flow  simultaneously  existing  subsonic  and  super¬ 
sonic  flows  along  with  accompanying  shock  waves  are  regarded  as 
characteristic.  From  the  viewpoint  of  mathematics,  the  properties  of 
transonic  flow  must  be  described  by  solving  "mixed"  differential 
equations  which  are  elliptical  in  the  subsonic  region  and  hyperboloid 
in  the  supersonic  region.  Equation  of  this  type  are  non-linear  and 
their  solution  normally  involve:.  ^continuous  surface  —  a  shock 
wave.  The  best  method  of  dealing  wi  '  transonic  flows  of  certain 
imbedding  shock  waves  is  the  finite  difference  relaxation  iteration 
method.  Although  the  first  to  use  the  relaxation  method  to  explain 
transonic  flows  was  Emmons^1],  but  since  he  used  the  Rankine-Hugoniot 
relationship  to  set  up  shock  waves,  therefore  widespread  application 
has  not  obtained.  The  finite  difference  relaxation  method  which 
automatically  calculates  the  shock  waves  which  was  put  forward  by 
Murman  and  Colet2]  has  opened  new  avenues  of  research  for  calculating 
transonic  flow  through  the  flow  field. 
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This  paper  is  based  on  the  characteristic  of  small  perturbations 
far  from  the  flow  about  a  body  (including  trailing  vortex)  to  rapidly 
attenuate,  causing  the  grid  to  tightly  flow  about  the  circumference 
of  the  object  without  affecting  the  accuracy  of  the  differential  oper¬ 
ator,  and  to  fully  utilize  accurate  boundary  conditions  of  the  distant 
flow  field,  particularly  using  the  hyperbolic  tangent  function  coord¬ 
inate  transformation  method,  an  infinitely  large  physical  space  is 
transformed  into  a  perfect  cube  with  side  lengths  of  [-1,  +1]  within 
which  the  finite  difference  calculations  are  carried  out. 

1.  TRANSONIC  SMALL  PERTURBATION  DIFFERENTIAL  EQUATIONS  AND 
BOUNDARY  CONDITIONS 

The  transonic  small  perturbation  differential  equation  disregard¬ 
ing  high-order  small  quantities  is: 

[l  -  ML  -  (r  +  0—  vw.Jvu...  +  -  o  (  1 ) 

for  convenience  of  calculation,  it  is  expressed  as  a  coordinate  tran¬ 
sformation,  let 

*  "  ■*■»/<■  —  yjb,*  —  S*MY»(r  + 

<p  -  !M2’(r  +  I)*/*.,:*1*]*, 

after  these  are  substituted  into  the  above  equation  and  arranged,  we 
obtain 

[X,  —  *  ,]*„  +  K&„  +  *„  —  o  ( 2  ) 


where 


x,  - 


1  -  ML 


K,  — 


(*Ml(r  +  a’faMLCr+OJ"* 

are  similar  parameters,  c  is  the  wing  root  half-chord,  b  is  the  half¬ 
span,  6  is  the  maximum  relative  thickness  of  the  airfoil. 


At  the  object  surface  with  —  /(*,?)  »  flow  condition 

~ •/*  is  satisfied,  and  a  is  the  angle  of  attack.  The 
Kutta  condition  is  satisfied  at  the  wing  trailing  edge.  Passing 
through  the  trailing  vortex,  and  are  both  continuous  but  is 
discontinuous,  its  jump  value  y,  ±o)  is  independent  of  x, 
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and  a*e  discontinuous.  At  the  boundary  of  the  distant  flow  field 
with  other  than  x  -+  00  ,  and  its  first  and  second  order  partial 
derivatives  of  x,  y,  and  2  are  all  equal  to  zero.  At  the  trailing 
vortex  with  x  ,  the  value  must  be  determined  on  the  basis  of 

r(y>  at  the  trailing  vortex.  Using  the  following  equation,  a  constant 
pressure  coefficient  can,  be  obtained  from  : 

79m 

e'~  M2’(r  +  i),A  **’  (3) 

and  after  integration,  the  various  aerodynamic  quantities  can  be 
obtained . 

2.  DIFFERENCE  SCHEME,  COORDINATE  TRANSFORMATION  AND  RELAXATION 
ITERATION 

Equation  (2)  is  suitable  for  finite  difference  iteration  solut¬ 
ions.  When  **<  ,  the  equation  becomes  elliptical  (subsonic  flow) 

and  the  central  difference  scheme  is  used;  when  9k  >  *4  .  the  equat¬ 
ion  become  hyperbolic  (supersonic  flow)  and  the  x  differential  uses  a 
one-sided  difference  operator.  Since  the  sonic  line  is  a  continuous 
transition  line,  no  special  provisions  must  be  made.  High  order  error 
terms  brought  on  by  replacement  of  differentials  by  difference  opera¬ 
tors,  give  rise  to  dissipative  effects;  consequently,  this  can  cause 
the  weak  shock  waves  whose  intensities  and  positions  were  previously 
unknown,  to  form  "naturally”  in  the  process  of  solution.  This  method 
is  called  the  shock  wave  capture  method.  In  order  to  allow  large 
changes  in  grid  density  in  the  flow  about  the  vicinity  of  the  object 
without  reducing  the  accuracy  of  the  solution,  we  employed  in  this 
paper  the  hyperbola  tangent  function  coordinate  transformation  method: 
$  —  thfcx,  n  —  th£f,  C  —  thf»  ,  and  the  infinitely  large  space  is  changed  into 
a  perfect  cube  with  side  lengths  Of  +1.  The  changes  are  a  one-to-one 
correspondence.  When  a,?,  f>\  (this  does  not  lose  its  generality), 
the  original  points  are  situated  in  the  flow  about  the  object  in  the 
center  of  the  calculated  space  and  occupy  76 %  of  the  space  of  the 
perfect  cube.  After  coordinate  conversion,  equation  (2)  can  be 
written  as: 

»,o  -  eot  k,  -  »o  -  r)*,no  -  -  2**1 

+  *J?0  -  lOIO  -  V)*,,  -  2^,1  +  f’(  1  -  C*)[(1  -  -  2Ce,l  -  0  <*»  > 
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Calculating  change  coefficient  linear  differential  equation  (U) 
is  not  nearly  so  difficult  as  equation  (2).  Here  we  omitted  solution 
of  approximate  potential  volume  integrals  for  finite  far-side  values, 
which  has  greatly  saved  computer  time.  This  partial  computation 
corresponds  to  the  time  required  to  solve  total  flow  field  iteration 
[33  • 

An  isometric  grid  is  marked  off  in  the  perfect  cube,  correspond¬ 
ing  to  a  non-isometric  grid  of  great  density  in  the  vicinity  of  the 
flow  about  the  object  in  the  physical  space.  For  each  grid  point  a 
difference  equation  is  established  and  the  linear  relaxation  iteration 
method  is  usedt1*3.  The  relaxation  lines  are  parallel  to  the  £  axis. 
The  potentials  for  points  in  the  proximity  of  the  relaxation  line  are 
regarded  as  known  values  for  solving  the  sets  of  difference  equations. 
After  the  potential  is  determined  for  the  nth  iteration,  ^<*+0 
calculated  on  the  basis  of  iteration  formula  +  O  — 

is  taken  as  the  initial  value,  replacing  ^H*on  the  relaxation  line. 

In  the  equation,  «a  is  the  relaxation  factor,  in  the  subsonic  region 
J,  £  m  <  Z  (super  relaxation),  and  in  the  supersonic  region  d<«o£j 
(low  relaxation).  As  soon  as  preliminary  estimates  of  the  ^  values 
are  made,  the  iteration  scan  can  be  carried  out.  Beginning  at  the 
section  with  greatest  value,  we  scan  from  the  upstream  boundary 
toward  the  downstream  boundary.  After  scanning  to  the  line  where 
i|*o  we  again  scan  one  section  at  a  time  along  the  i\  value  or  small 
direction  until  the  entire  flow  field  has  been  scanned^then  we  con¬ 
sider  that  one  iteration  is  complete.  After  each  iteration,  we  use 
the  values  for  this  iteration  to  determine  the  circulation  distribu¬ 
tion.  The  ^  obtained  at  the  Trefftz  plane  is  taken  as  the  boundary 
value  of  the  lower  iteration.  After  iteration  is  carried  out  a  cer¬ 
tain  number  of  times,  if  every  point  in  the  flow  field  satisfies 

< •  (c  is  control  precision,  generally  taken  as  10~5)t  then 
we  consider  that  the  solution  is  convergent.  The  optimum  value  range 
of  6»  confirmed  by  many  numerical  tests:  super  relaxation  factor  1.3 
&*»£!» f*  for  example,  w  =  1.7  saves  1/3  the  computer  time  of«v  s  1; 
the  low  relaxation  factor  0.5Sc*»£».i.  In  the  process  of  the  calcul¬ 
ations,  if  all  the  points  on  the  entire  relaxation  line  are  subsonic 


points,  then  the  set  of  difference  equations  will  be  linear;  if 
supersonic  points  occur  on  the  relaxation  line,  then  they  will  be 
nonlinear.  Now,  relatively  close  values  should  be  used  for  solving 
initial  value  iteration.  For  example,  results  obtained  when  we  begin 
with  a  certain  subcritical  Mae  are  used  to  calculate  initial  values 
for  relatively  high  M  «e  (the  Mach  number  climbing  method).  Numerical 
tests  indicate  that  so  long  as  the  <*>  value  initially  selected  is  low, 
we  can  still  start  with  Mo*  =  0  to  calculate  ^  under  high  subcritical 
M*».  Under  supercritical  conditions,  if  the  Mach  number  climb  inter¬ 
val  is  too  great,  such  as  a  Maio.j,  there  is  the  possibility  of 
resulting  in  iteration  divergence.  In  this  case  the  climb  interval 
should  immediately  be  decresed  and  iteration  continued.  Computer 
time  is  directly  proportional  to  the  (1  +  r)th  power  of  the  grid 
point  o.S  and  directly  proportional  to  the  logarithm  of  con¬ 

trol  precision  £,  . 

3.  CALCULATED  RESULTS 

In  this  paper  we  have  calculated  pressure  distributions  for 
rectangular  wings  with  a  6%  double  circular-arc  airfoil  and  an  aspect 
ratio  of  4,  wings  with  a  constant  chord  and  30°  sweep-back,  as  well 
as  rectangular  wings  with  a  5%  double  circular-are  airfoil  and  an 
aspect  ratio  of  3.  From  the  illustrations  it  can  be  seen  that  our 
calculated  results  coincide  well  with  the  test  results  in  [5]  and 
[8]  and  coincide  with  the  results  calculated  in  [6]  and  [7].  The 
results  of  subcritical  to  supercritical  and  zero  angle  of  attack  to 
having  an  angle  of  attack  coincide  very  well  with  test  results.  This 
is  because  of  difference  scheme  of  low  precision  was  used  in  the 
supersonic  region  and  the  flow  in  the  vicinity  of  the  leading  edge 
and  in  the  shock  wave  region  undergoes  rather  drastic  changes.  The 
method  employed  in  this  paper  decreases  computer  time  by  approximate¬ 
ly  one-half  that  of  [6].  Taking  the  results  in  Fig.  1  as  an  example, 
the  19X10X19  grid  points  used  in  this  paper  involved  a  computer  time 
of  18  minutes  on  a  655  machine,  but  the  50X50X20  grid  points  used  in 
[6]  involved  a  computer  time  of  30  minutes  on  an  IBM  360/67.  Surface 
relaxation  was  also  tested  in  this  paper.  Surface  relaxation  was 
faster  than  linear  relaxation. 


The  difficulty  encountered  in  calculating  three-dimensional 
transonic  flow  is  how  to  deal  with  boundary  conditions  for  complicat¬ 
ed  profiles.  Of  particular  difficulty  are  linearizing  boundary  cond¬ 
itions  for  working  out  solutions  in  the  vicinity  of  a  blunt  leading 
edge  and  improving  calculations  for  swept-back  shock  waves. 


Fig.  1.  Constant  pressure  distribution 
a  6%  double  circular-arc  airfoil  and  an 

O  test  values  t5]f  . . —  calculated 

the  half-span  section,  M  oe  =  0.806,  <t 


for  a  rectangular  wing  with 
aspect  ratio  of  4. 

values  for  this  paper,  y  is 
=  0°,  subcritical 


Fig.  2.  Constant  pressure  distrib¬ 
ution  for  a  rectangular  wing  with  a 
6 %  double  circular-arc  profile  and 
an  aspect  ratio  of  4. 


o  test  values^],  y  =  0.  A  test 
valueu  (two-dimensional ) ] t 

-  calculated  value  for  this 

paper,  M  «•  =  0.908,  A:  0°,  super¬ 
critical  . 


Fig.  3-  Constant  pressure  dist¬ 
ribution  for  a  rectangular  wing 
a  6 %  double  circular-arc  airfoil 
and  an  aspect  ratio  of  4. 

®  upperTgurface  tes£  values 1 5 3 
a  lowerj 

—  calculated  values  for  this 
paper,  M  —  s  0.908,  «  =  1°, 
supercritical. 
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Fig.  Wing  root  cross-sectional  pressure  distribution  for  a  30° 
swept-back  wing  with  a  6%  double  circular-arc  profile  and  an  aspect 
ratio  of  4. 
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Fig.  5.  Constant  pressure  distribution  for  a  rectangular  wing  with 
a  5%  double  circular-arc  airfoil  and  an  aspect  ratio  of  3* 

o  upper  surface  test  values^®},  v  lower  surface  test  values, 

- - — .  calculated  values  for  this  paper 

Moo  =  0.9,  «  =  5°,  y  =  0  supercritical 
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